// Стандартные
#define _USE_MATH_DEFINES   // Включения определений мат. констант (например числа пи), не прописаных в стандарте C++.
#include <cmath>
// QtCore
#include <QtCore/qmath.h>
#include <QtCore/QDir>
#include <QtCore/QTextCodec>
#include <QtCore/QDataStream>
#include <QtCore/QStringBuilder>
// QtGui
#include <QtGui/QApplication>
// Локальные
#include "settings.h"
#include "linal/los.h"
#include "linal/bcg.h"
#include "abstractsolver.h"
#include "fem/trianglefem.h"
#include "approximation/interpolation.h"

//================================================================================================================================
// Значение производной по направлению в точке (узле конечноэлементной сетки).
//================================================================================================================================
double AbstractSolver::derivation(const linal::Vector &f, int point, DiffVar var, DiffMethod method)
{
    double dx;
    geometry::_2d::Point p;

    if(method == AbstractSolver::LDiff)
        dx = - m_DiffStep;
    else
        dx = m_DiffStep;

    if(var == AbstractSolver::R)
        p = m_Net.point(point).shifted(dx, 0);
    else
        p = m_Net.point(point).shifted(0, dx);

    int triangle = m_NetSearcher->findTriangle(m_Net, p);
    double df = fem::triangles::functionValue(triangle, m_Net, f, p) - f[point];

    return df/dx;
}

//================================================================================================================================
//================================================================================================================================
AbstractSolver::AbstractSolver() : m_NetSearcher(NULL)
{
    //****************************************************************************************************************************
    // Создание и очистка директории для вывода результатов.
    //****************************************************************************************************************************
    m_OutputDir = Settings::instance()->setting("OutputPath").toString();
    if(!QDir::isAbsolutePath(m_OutputDir))
        m_OutputDir = QDir::cleanPath(m_OutputDir).prepend(qApp->applicationDirPath() + "/");
    m_InputDir = Settings::instance()->setting("InputPath").toString();
    if(!QDir::isAbsolutePath(m_InputDir))
        m_InputDir = QDir::cleanPath(m_InputDir).prepend(qApp->applicationDirPath() + "/");

    if(!QDir(m_OutputDir).exists())
        QDir(qApp->applicationDirPath()).mkpath(m_OutputDir);
    if(!QDir(m_InputDir).exists())
        QDir(qApp->applicationDirPath()).mkpath(m_InputDir);
    // TODO: Очистка директории.

    //****************************************************************************************************************************
    // Считывание параметров временной сетки.
    //****************************************************************************************************************************
    m_NumOfTimeLayers = Settings::instance()->setting("Time/TimeLayers").toInt();
    m_TimeMoment = Settings::instance()->setting("Time/TimeMoment").toInt();
    m_TimeStep = Settings::instance()->setting("Time/TimeStep").toDouble();

    if(Settings::instance()->setting("Time/Stationarity").toBool())
        m_SubTasks |= AbstractSolver::Stationarity;

    //****************************************************************************************************************************
    // Считывание физических и геометрических параметров задачи. Перевод скорости из см/ч в м/с.
    //****************************************************************************************************************************
    m_Radius = Settings::instance()->setting("Physics/CrucibleRadius").toDouble();
    m_TMax = Settings::instance()->setting("Physics/MaxTemperature").toDouble();
    m_TRod = Settings::instance()->setting("Physics/RodTemperature").toDouble();
    m_TGrad = Settings::instance()->setting("Physics/TemperatureGradient").toDouble();
    m_TGrad *= 100.;    //К/см -> К/м
    m_VRot = Settings::instance()->setting("Physics/AngularSpeed").toDouble();
    m_V = Settings::instance()->setting("Physics/CrucibleSpeed").toDouble();
    m_V /= 100.;        //см/ч -> м/ч
    m_V /= 3600.;       //м/ч -> м/с

    m_MinCrystalHeight = Settings::instance()->setting("Math/MinCrystalHeight").toDouble();
    //****************************************************************************************************************************
    // Считывание параметров решателей.
    //****************************************************************************************************************************
    m_MaxiterSLAE = Settings::instance()->setting("Math/SlaeIters").toInt();
    m_MaxiterNonlin = Settings::instance()->setting("Math/NonlinIters").toInt();

    m_EpsSLAE = Settings::instance()->setting("Math/SlaeResidual").toDouble();
    m_EpsNonlin = Settings::instance()->setting("Math/NonlinResidual").toDouble();

    //****************************************************************************************************************************
    // Инициализация решателей.
    //****************************************************************************************************************************
    m_Solvers.append(new linal::LOS(m_MaxiterSLAE, m_EpsSLAE));
    m_Solvers.append(new linal::BCG(m_MaxiterSLAE, m_EpsSLAE));

    //****************************************************************************************************************************
    // Считывание полуширины фронта.
    //****************************************************************************************************************************
    m_FrontHalfWidth = Settings::instance()->setting("Math/FrontHalfwidth").toDouble();

    //****************************************************************************************************************************
    // Считывание шага производной.
    //****************************************************************************************************************************
    m_DiffStep = Settings::instance()->setting("Math/DiffStep").toDouble();

    //****************************************************************************************************************************
    // Считывание параметра релаксации.
    //****************************************************************************************************************************
    m_W = Settings::instance()->setting("Math/RelaxFactor").toDouble();

    //****************************************************************************************************************************
    // Считывание сервисных настроек программы.
    //****************************************************************************************************************************
    m_OutputStep = Settings::instance()->setting("Solver/OtputStep").toInt();
    m_ResumeFrom = Settings::instance()->setting("Solver/ResumeFrom").toInt();

    //****************************************************************************************************************************
    // Считывание флагов.
    //****************************************************************************************************************************
    if(Settings::instance()->setting("Solver/Solidification").toBool())
        m_SubTasks |= AbstractSolver::Solidification;
    if(Settings::instance()->setting("Solver/Convection").toBool())
        m_SubTasks |= AbstractSolver::Convection;
    if(Settings::instance()->setting("Solver/Rotation").toBool())
        m_SubTasks |= AbstractSolver::Rotation;
    if(Settings::instance()->setting("Solver/Concentration").toBool())
        m_SubTasks |= AbstractSolver::Concentration;
    if(Settings::instance()->setting("Solver/Debug").toBool())
        m_SubTasks |= AbstractSolver::Debug;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::bin2txt()
{
    //****************************************************************************************************************************
    // Чтение и обработка настроек.
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_dir(m_InputDir);
    
    QFile in_file;
    QDataStream istream(&in_file);
    istream.setVersion(QDataStream::Qt_4_7);

    // Файл для записи результата.
    QFile out_file;
    QTextStream ostream(&out_file);
    ostream.setRealNumberPrecision(12);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Bin2Txt") << endl;
    foreach(QString filename, in_dir.entryList(QStringList("*nvtr.dat"), QDir::Files | QDir::CaseSensitive))
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();
        in_file.setFileName(m_InputDir + "/" + filename);
        out_file.setFileName(m_OutputDir + "/" + filename);
        in_file.open(QFile::ReadOnly);
        out_file.open(QFile::WriteOnly | QFile::Text);
        while(!istream.atEnd())
        {
            int points[3], material;
            istream >> points[0] >> points[1] >> points[2] >> material;
            ostream << points[0] << " " << points[1] << " " << points[2] << " " << material << endl;
        }
        in_file.close();
        out_file.close();

        QString filename_rz = filename.left(filename.size() - 8) + "rz.dat";
        cout << QObject::tr("Файл") << " " << filename_rz << "\r";
        cout.flush();
        in_file.setFileName(m_InputDir + "/" + filename_rz);
        out_file.setFileName(m_OutputDir + "/" + filename_rz);
        in_file.open(QFile::ReadOnly);
        out_file.open(QFile::WriteOnly | QFile::Text);
        while(!istream.atEnd())
        {
            double r, z;
            istream >> r >> z;
            ostream << r << " " << z << endl;
        }
        in_file.close();
        out_file.close();

        QString filename_t = filename.left(filename.size() - 8) + "T.dat";
        cout << QObject::tr("Файл") << " " << filename_t << "\r";
        cout.flush();
        in_file.setFileName(m_InputDir + "/" + filename_t);
        out_file.setFileName(m_OutputDir + "/" + filename_t);
        in_file.open(QFile::ReadOnly);
        out_file.open(QFile::WriteOnly | QFile::Text);
        while(!istream.atEnd())
        {
            double r, z, f;
            istream >> r >> z >> f;
            ostream << r << " " << z << " " << f << endl;
        }
        in_file.close();
        out_file.close();

        if(m_SubTasks.testFlag(AbstractSolver::Concentration))
        {
            QString filename_c = filename.left(filename.size() - 8) + "C.dat";
            cout << QObject::tr("Файл") << " " << filename_c << "\r";
            cout.flush();
            in_file.setFileName(m_InputDir + "/" + filename_c);
            out_file.setFileName(m_OutputDir + "/" + filename_c);
            in_file.open(QFile::ReadOnly);
            out_file.open(QFile::WriteOnly | QFile::Text);
            while(!istream.atEnd())
            {
                double r, z, f;
                istream >> r >> z >> f;
                ostream << r << " " << z << " " << f << endl;
            }
            in_file.close();
            out_file.close();
        }

        if(m_SubTasks.testFlag(AbstractSolver::Convection))
        {
            QString filename_vr = filename.left(filename.size() - 8) + "Vr.dat";
            cout << QObject::tr("Файл") << " " << filename_vr << "\r";
            cout.flush();
            in_file.setFileName(m_InputDir + "/" + filename_vr);
            out_file.setFileName(m_OutputDir + "/" + filename_vr);
            in_file.open(QFile::ReadOnly);
            out_file.open(QFile::WriteOnly | QFile::Text);
            while(!istream.atEnd())
            {
                double r, z, f;
                istream >> r >> z >> f;
                ostream << r << " " << z << " " << f << endl;
            }
            in_file.close();
            out_file.close();


            QString filename_vz = filename.left(filename.size() - 8) + "Vz.dat";
            cout << QObject::tr("Файл") << " " << filename_vz << "\r";
            cout.flush();
            in_file.setFileName(m_InputDir + "/" + filename_vz);
            out_file.setFileName(m_OutputDir + "/" + filename_vz);
            in_file.open(QFile::ReadOnly);
            out_file.open(QFile::WriteOnly | QFile::Text);
            while(!istream.atEnd())
            {
                double r, z, f;
                istream >> r >> z >> f;
                ostream << r << " " << z << " " << f << endl;
            }
            in_file.close();
            out_file.close();

            if(m_SubTasks.testFlag(AbstractSolver::Rotation))
            {
                QString filename_vf = filename.left(filename.size() - 8) + "Vf.dat";
                cout << QObject::tr("Файл") << " " << filename_vf << "\r";
                cout.flush();
                in_file.setFileName(m_InputDir + "/" + filename_vf);
                out_file.setFileName(m_OutputDir + "/" + filename_vf);
                in_file.open(QFile::ReadOnly);
                out_file.open(QFile::WriteOnly | QFile::Text);
                while(!istream.atEnd())
                {
                    double r, z, f;
                    istream >> r >> z >> f;
                    ostream << r << " " << z << " " << f << endl;
                }
                in_file.close();
                out_file.close();
            }

            QString filename_psi = filename.left(filename.size() - 8) + "Psi.dat";
            cout << QObject::tr("Файл") << " " << filename_psi << "\r";
            cout.flush();
            in_file.setFileName(m_InputDir + "/" + filename_psi);
            out_file.setFileName(m_OutputDir + "/" + filename_psi);
            in_file.open(QFile::ReadOnly);
            out_file.open(QFile::WriteOnly | QFile::Text);
            while(!istream.atEnd())
            {
                double r, z, f;
                istream >> r >> z >> f;
                ostream << r << " " << z << " " << f << endl;
            }
            in_file.close();
            out_file.close();

            QString filename_omega = filename.left(filename.size() - 8) + "Omega.dat";
            cout << QObject::tr("Файл") << " " << filename_omega << "\r";
            cout.flush();
            in_file.setFileName(m_InputDir + "/" + filename_omega);
            out_file.setFileName(m_OutputDir + "/" + filename_omega);
            in_file.open(QFile::ReadOnly);
            out_file.open(QFile::WriteOnly | QFile::Text);
            while(!istream.atEnd())
            {
                double r, z, f;
                istream >> r >> z >> f;
                ostream << r << " " << z << " " << f << endl;
            }
            in_file.close();
            out_file.close();
        }
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::postproc()
{
    if(Settings::instance()->setting("Postproc/Bin2Txt").toBool())
        bin2txt();
    if(Settings::instance()->setting("Postproc/SpeedProfiles").toBool())
        speedProfiles();
    if(Settings::instance()->setting("Postproc/TemperatureProfiles").toBool())
        temperatureProfiles();
    if(Settings::instance()->setting("Postproc/CrystalMass").toBool())
        crystalMass();
    if(Settings::instance()->setting("Postproc/BuildFront").toBool())
        crystallizationFront();
    if(Settings::instance()->setting("Postproc/KineticEnergy").toBool())
        kineticEnergy();
    if(Settings::instance()->setting("Postproc/LocalConductiveHeatFlow").toBool())
        localConductiveHeatFlow();
    if(Settings::instance()->setting("Postproc/LocalConvectiveHeatFlow").toBool())
        localConvectiveHeatFlow();
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::crystalMass()
{
    //****************************************************************************************************************************
    // Чтение и обработка настроек.
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_dir(m_InputDir);

    // Файл для записи результата.
    QFile out_file(m_OutputDir + "/crystal mass.dat");
    QTextStream ostream(&out_file);

    if(out_file.exists())
        out_file.remove();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Масса кристалла") << endl;
    foreach(QString filename, in_dir.entryList(QStringList("*nvtr.dat"), QDir::Files | QDir::CaseSensitive))
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        QString rz_filename = in_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "rz.dat";
        QString nvtr_filename = in_dir.absolutePath() + "/" + filename;
        m_Net.readFromFiles(rz_filename, nvtr_filename, fem::TriangleNet::NET_T);

        //************************************************************************************************************************
        // По информации о сетке находим объем кристалла, а затем его массу.
        //************************************************************************************************************************
        double mass = 0;
        for(int n = 0; n < m_Net.numOfElements(); n++)
        {
            if(m_Materials[m_Net.material(n)].state() == Material::Liquid)
            {
                double detD = (m_Net.point(n,1).x() - m_Net.point(n,0).x())*(m_Net.point(n,2).y() - m_Net.point(n,0).y()) -
                        (m_Net.point(n,2).x() - m_Net.point(n,0).x())*(m_Net.point(n,1).y() - m_Net.point(n,0).y());

                for(int i = 0; i < 3; i++)
                    mass += 2*M_PI*m_Materials[m_Net.material(n)].property(Material::Density)*m_Net.point(n,i).x()*qAbs(detD)/6.;
            }
        }

        //************************************************************************************************************************
        // Записываем результат в файл.
        //************************************************************************************************************************
        out_file.open(QFile::Append | QFile::Text);
        ostream << filename.remove("nvtr.dat").toInt()*m_TimeStep/60. << " " << mass << endl;
        out_file.close();
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::kineticEnergy()
{
    //****************************************************************************************************************************
    // Чтение и обработка настроек.
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_dir(m_InputDir);

    // Файл для записи результата.
    QFile out_file1(m_OutputDir + "/kinetic energy.dat");
    QFile out_file2(m_OutputDir + "/specific kinetic energy.dat");
    QTextStream ostream1(&out_file1);
    QTextStream ostream2(&out_file2);

    if(out_file1.exists())
        out_file1.remove();
    if(out_file2.exists())
        out_file2.remove();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Кинетическая энергия расплава") << endl;
    foreach(QString filename, in_dir.entryList(QStringList("*nvtr.dat"), QDir::Files | QDir::CaseSensitive))
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        QString rz_filename = in_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "rz.dat";
        QString nvtr_filename = in_dir.absolutePath() + "/" + filename;
        m_Net.readFromFiles(rz_filename, nvtr_filename, fem::TriangleNet::NET_T);

        m_Vr.readFromFile(in_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "Vr.dat", linal::Vector::FEM_T);
        m_Vz.readFromFile(in_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "Vz.dat", linal::Vector::FEM_T);

        linal::Vector V(m_Net.numOfPoints());
        for(int i = 0; i < V.dimension(); i++)
            V[i] = m_Vr[i]*m_Vr[i] + m_Vz[i]*m_Vz[i];

        //************************************************************************************************************************
        //************************************************************************************************************************
        double energy = 0;
        double mass = 0;
        double C[3][3] = {{1./12., 1./24., 1./24.}, {1./24., 1./12., 1./24.}, {1./24., 1./24., 1./12.}};

        for(int n = 0; n < m_Net.numOfElements(); n++)
        {
            if(m_Materials[m_Net.material(n)].state() == Material::Liquid)
            {
                double detD = (m_Net.point(n,1).x() - m_Net.point(n,0).x())*(m_Net.point(n,2).y() - m_Net.point(n,0).y()) -
                        (m_Net.point(n,2).x() - m_Net.point(n,0).x())*(m_Net.point(n,1).y() - m_Net.point(n,0).y());
                double density = m_Materials[m_Net.material(n)].property(Material::Density);

                for(int i = 0; i < 3; i++)
                    for(int j = 0; j < 3; j++)
                        energy += 2*M_PI*density*m_Net.point(n,i).x()*(V[m_Net.pointIndex(n, j)]*C[i][j])*qAbs(detD)/2.;

                for(int i = 0; i < 3; i++)
                    mass += 2*M_PI*density*m_Net.point(n,i).x()*qAbs(detD)/6.;
            }
        }

        //************************************************************************************************************************
        // Записываем результат в файл.
        //************************************************************************************************************************
        out_file1.open(QFile::Append | QFile::Text);
        ostream1 << filename.remove("nvtr.dat").toInt()*m_TimeStep/60. << " " << energy << endl;
        out_file1.close();

        out_file2.open(QFile::Append | QFile::Text);
        ostream2 << filename.remove("nvtr.dat").toInt()*m_TimeStep/60. << " " << energy/mass << endl;
        out_file2.close();
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::crystallizationFront()
{
    //****************************************************************************************************************************
    // Чтение и обработка настроек.
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_dir(m_InputDir);

    // Директория для записи результата.
    QDir out_dir(m_OutputDir + "/Front");
    if(!out_dir.exists())
        out_dir.mkpath(out_dir.absolutePath());

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Фронт кристаллизации") << endl;
    foreach(QString filename, in_dir.entryList(QStringList("*nvtr.dat"), QDir::Files | QDir::CaseSensitive))
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        QString rz_filename = in_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "rz.dat";
        QString nvtr_filename = in_dir.absolutePath() + "/" + filename;
        m_Net.readFromFiles(rz_filename, nvtr_filename, fem::TriangleNet::NET_T);

        m_T.readFromFile(in_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "T.dat", linal::Vector::FEM_T);

        //************************************************************************************************************************
        //Выделяем фронт кристаллизации.
        //************************************************************************************************************************
        QList<geometry::_2d::Point> front = identifyFrontExt(1685 + m_FrontHalfWidth);

        //************************************************************************************************************************
        //************************************************************************************************************************
        QFile file(out_dir.absolutePath() + "/" + filename.left(filename.size() - 8) + "Front.dat");
        QTextStream ostream(&file);

        file.open(QFile::WriteOnly | QFile::Text);
        foreach(geometry::_2d::Point point, front)
            ostream << point << endl;
        file.close();
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::localConvectiveHeatFlow()
{
    // Директория с входными данными.
    QDir in_path(m_InputDir);
    // Маска файла/файлов для обработки.
    QStringList filename_mask(Settings::instance()->setting("LocalConvectiveHeatFlow/FileMask").toString() + "T.dat");
    // Директория для записи результата.
    QDir out_path(m_OutputDir);

    // Получаем список файлов для обработки.
    QStringList filenames = in_path.entryList(filename_mask, QDir::CaseSensitive | QDir::Files);

    //Создаем структуру директорий для записи результата.
    QDir qz_r_out_path(out_path.absolutePath() + "/Qz(r) Conv");
    QDir qz_z_out_path(out_path.absolutePath() + "/Qz(z) Conv");
    QDir qr_z_out_path(out_path.absolutePath() + "/Qr(z) Conv");
    QDir qz_r_s_out_path(out_path.absolutePath() + "/Qz(r) Conv/front");

    if(!qz_r_out_path.exists())
        qz_r_out_path.mkpath(qz_r_out_path.absolutePath());

    if(!qz_z_out_path.exists())
        qz_z_out_path.mkpath(qz_z_out_path.absolutePath());

    if(!qr_z_out_path.exists())
        qr_z_out_path.mkpath(qr_z_out_path.absolutePath());

    if(!qz_r_s_out_path.exists())
        qz_r_s_out_path.mkpath(qz_r_s_out_path.absolutePath());

    //****************************************************************************************************************************
    // Получение задаваемых извне координат сечений.
    //****************************************************************************************************************************
    QVariantList r_sections_temp = Settings::instance()->setting("LocalConvectiveHeatFlow/R").toList();
    QVariantList z_sections_temp = Settings::instance()->setting("LocalConvectiveHeatFlow/Z").toList();

    QList<double> r_sections;
    QList<double> z_sections;

    for(int i = 0; i < r_sections_temp.size(); i++)
        r_sections.append(r_sections_temp[i].toDouble());

    for(int i = 0; i < z_sections_temp.size(); i++)
        z_sections.append(z_sections_temp[i].toDouble());

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    double density = m_Materials[1].property(Material::Density);
    double heat_capacity = m_Materials[1].property(Material::HeatCapacity);

    //****************************************************************************************************************************
    // Обработка файлов. Генерация тепловых потоков.
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Локальный конвективный тепловой поток") << endl;
    foreach(QString filename, filenames)
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        //************************************************************************************************************************
        //************************************************************************************************************************
        m_T.readFromFile(in_path.absolutePath() + "/" + filename, linal::Vector::FEM_T);
        m_Vr.readFromFile(in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "Vr.dat", linal::Vector::FEM_T);
        m_Vz.readFromFile(in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "Vz.dat", linal::Vector::FEM_T);
        m_Net.readFromFiles(in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "rz.dat",
                            in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "nvtr.dat", fem::TriangleNet::NET_T);

        m_NetSearcher->analyzeNet(m_Net);

        //************************************************************************************************************************
        //************************************************************************************************************************
        foreach(double r, r_sections)
        {
            if(!qz_z_out_path.exists(QString::number(r)))
                qz_z_out_path.mkpath(QString::number(r));

            if(!qr_z_out_path.exists(QString::number(r)))
                qr_z_out_path.mkpath(QString::number(r));

            int n;
            double z_min, z_max;
            limitsForZ(r, &z_min, &z_max, &n, true);

            QFile file1(qz_z_out_path.absolutePath() + "/" + QString::number(r) + "/" + filename);
            QFile file2(qr_z_out_path.absolutePath() + "/" + QString::number(r) + "/" + filename);
            QTextStream ostream1(&file1);
            QTextStream ostream2(&file2);

            file1.open(QFile::WriteOnly | QFile::Text);
            file2.open(QFile::WriteOnly | QFile::Text);
            for(int j = 0; j < n; j++)
            {
                double z = z_min + j*(z_max - z_min)/(n - 1);
                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));

                double t = fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(r, z));
                double vz = fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(r, z));
                double vr = fem::triangles::functionValue(triangle, m_Net, m_Vr, geometry::_2d::Point(r, z));

                ostream1 << density*heat_capacity*vz*t << " " << z << endl;
                ostream2 << density*heat_capacity*vr*t << " " << z << endl;
            }
            file1.close();
            file2.close();
        }

        //************************************************************************************************************************
        //************************************************************************************************************************
        foreach(double z, z_sections)
        {
            if(!qz_r_out_path.exists(QString::number(z)))
                qz_r_out_path.mkpath(QString::number(z));

            int n;
            double r_min, r_max;
            limitsForR(z, &r_min, &r_max, &n, true);

            QFile file(qz_r_out_path.absolutePath() + "/" + QString::number(z) + "/" + filename);
            QTextStream ostream(&file);

            file.open(QFile::WriteOnly | QFile::Text);
            for(int i = 0; i < n; i++)
            {
                double r = r_min + i*(r_max - r_min)/(n - 1);
                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));

                double t = fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(r, z));
                double vz = fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(r, z));

                ostream << r << " " << density*heat_capacity*vz*t << endl;
            }
            file.close();
        }

        //************************************************************************************************************************
        //Распределение осевого потока тепла в радиальном направлении (прямая над фронтом).
        //************************************************************************************************************************
        if(Settings::instance()->setting("LocalConvectiveHeatFlow/FlowAboveFront").toBool())
        {
            QList<geometry::_2d::Point> front_plus = identifyFrontExt(1685 + m_FrontHalfWidth);

            //********************************************************************************************************************
            // Выводим результат.
            //********************************************************************************************************************
            QFile file(qz_r_s_out_path.absolutePath() + "/" + filename);
            QTextStream ostream(&file);
            file.open(QFile::WriteOnly | QFile::Text);

            for(int i = 0; i < front_plus.size(); i++)
            {
                double r(front_plus[i].x());
                double z(front_plus[Settings::instance()->setting("LocalConvectiveHeatFlow/CurvedSectionAboveFront").toBool()? i : 0].y() + 0.00125);

                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));

                double t = fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(r, z));
                double vz = fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(r, z));

                ostream << r << " " << density*heat_capacity*vz*t << endl;
            }
            file.close();
        }
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::localConductiveHeatFlow()
{
    //****************************************************************************************************************************
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_path(m_InputDir);
    // Маска файла/файлов для обработки.
    QStringList filename_mask(Settings::instance()->setting("LocalConductiveHeatFlow/FileMask").toString() + "T.dat");
    // Директория для записи результата.
    QDir out_path(m_OutputDir);

    // Получаем список файлов для обработки.
    QStringList filenames = in_path.entryList(filename_mask, QDir::CaseSensitive | QDir::Files);

    //Создаем структуру директорий для записи результата.
    QDir qz_z_out_path(out_path.absolutePath() + "/Qz(z) Cond");
    QDir qz_r_out_path(out_path.absolutePath() + "/Qz(r) Cond");
    QDir qr_z_w_out_path(out_path.absolutePath() + "/Qr(z) Cond/wall");
    QDir qz_r_s_out_path(out_path.absolutePath() + "/Qz(r) Cond/front");

    if(!qz_z_out_path.exists())
        qz_z_out_path.mkpath(qz_z_out_path.absolutePath());

    if(!qz_r_out_path.exists())
        qz_r_out_path.mkpath(qz_r_out_path.absolutePath());

    if(!qr_z_w_out_path.exists())
        qr_z_w_out_path.mkpath(qr_z_w_out_path.absolutePath());

    if(!qz_r_s_out_path.exists())
        qz_r_s_out_path.mkpath(qz_r_s_out_path.absolutePath());

    //****************************************************************************************************************************
    // Получение задаваемых извне координат сечений.
    //****************************************************************************************************************************
    QVariantList r_sections_temp = Settings::instance()->setting("LocalConductiveHeatFlow/R").toList();
    QVariantList z_sections_temp = Settings::instance()->setting("LocalConductiveHeatFlow/Z").toList();

    QList<double> r_sections;
    QList<double> z_sections;

    for(int i = 0; i < r_sections_temp.size(); i++)
        r_sections.append(r_sections_temp[i].toDouble());

    for(int i = 0; i < z_sections_temp.size(); i++)
        z_sections.append(z_sections_temp[i].toDouble());

    //****************************************************************************************************************************
    // Обработка файлов. Генерация тепловых потоков.
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Локальный кондуктивный тепловой поток") << endl;
    foreach(QString filename, filenames)
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        //************************************************************************************************************************
        //************************************************************************************************************************
        m_T.readFromFile(in_path.absolutePath() + "/" + filename, linal::Vector::FEM_T);
        m_Net.readFromFiles(in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "rz.dat",
                            in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "nvtr.dat", fem::TriangleNet::NET_T);

        //************************************************************************************************************************
        //Анализируем сетку и составляем таблицу, сокращающую поиск по сетке.
        //************************************************************************************************************************
        m_NetSearcher->analyzeNet(m_Net);

        //************************************************************************************************************************
        //Выделяем фронт кристаллизации.
        //************************************************************************************************************************
        QList<geometry::_2d::Point> front_plus = identifyFrontExt(1685 + m_FrontHalfWidth);

        //************************************************************************************************************************
        //Распределение осевого потока тепла по высоте тигля (расплав).
        //************************************************************************************************************************
        foreach(double r, r_sections)
        {
            //********************************************************************************************************************
            //********************************************************************************************************************
            QList<int> triangles;
            QList<double> temperatures;

            int n;
            double z_min, z_max;
            limitsForZ(r, &z_min, &z_max, &n, false);

            for(int j = 0; z_min + j*m_DiffStep < z_max; j++)
            {
                double z = z_min + j*m_DiffStep;
                triangles.append(m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z)));
                temperatures.append(fem::triangles::functionValue(triangles.last(), m_Net, m_T, geometry::_2d::Point(r, z)));
            }

            //********************************************************************************************************************
            //********************************************************************************************************************
            if(!qz_z_out_path.exists(QString::number(r)))
                qz_z_out_path.mkdir(QString::number(r));
            QFile file(qz_z_out_path.absolutePath() + "/" + QString::number(r) + "/" + filename);
            QTextStream ostream(&file);
            file.open(QFile::WriteOnly | QFile::Text);

            for(int i = 0; i < temperatures.size() - 1; i++)
            {
                if(m_Net.material(triangles[i]) == m_Net.material(triangles[i + 1]))
                {
                    int material = m_Net.material(triangles[i]);
                    ostream << m_Materials[material].property(Material::HeatConductivity)*(temperatures[i + 1] - temperatures[i])/m_DiffStep;
                    ostream << " " << z_min + i*m_DiffStep << endl;
                }
            }
        }

        //************************************************************************************************************************
        //Распределение осевого потока тепла в радиальном направлении (расплав).
        //************************************************************************************************************************
        foreach(double z, z_sections)
        {
            //********************************************************************************************************************
            //********************************************************************************************************************
            QList<int> triangles1, triangles2;
            QList<double> temperatures1, temperatures2;

            int n;
            double r_min, r_max;
            limitsForR(z, &r_min, &r_max, &n, false);

            for(int i = 0; r_min + i*m_DiffStep < r_max; i++)
            {
                double r = r_min + i*m_DiffStep;
                triangles1.append(m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z)));
                temperatures1.append(fem::triangles::functionValue(triangles1.last(), m_Net, m_T, geometry::_2d::Point(r, z)));

                triangles2.append(m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z + m_DiffStep)));
                temperatures2.append(fem::triangles::functionValue(triangles2.last(), m_Net, m_T, geometry::_2d::Point(r, z + m_DiffStep)));
            }

            //********************************************************************************************************************
            //********************************************************************************************************************
            if(!qz_r_out_path.exists(QString::number(z)))
                qz_r_out_path.mkdir(QString::number(z));
            QFile file(qz_r_out_path.absolutePath() + "/" + QString::number(z) + "/" + filename);
            QTextStream ostream(&file);
            file.open(QFile::WriteOnly | QFile::Text);

            for(int i = 0; i < temperatures1.size(); i++)
            {
                if(m_Net.material(triangles1[i]) == m_Net.material(triangles2[i]))
                {
                    int material = m_Net.material(triangles1[i]);
                    ostream << r_min + i*m_DiffStep << " ";
                    ostream << m_Materials[material].property(Material::HeatConductivity)*(temperatures2[i] - temperatures1[i])/m_DiffStep  << endl;
                }
            }
        }

        //************************************************************************************************************************
        //Распределение радиального потока тепла по высоте тигля (боковая стенка).
        //************************************************************************************************************************
        if(Settings::instance()->setting("LocalConductiveHeatFlow/FlowFromWall").toBool())
        {
            QFile file(qr_z_w_out_path.absolutePath() + "/" + filename);
            QTextStream ostream(&file);
            file.open(QFile::WriteOnly | QFile::Text);

            int n;
            double z_min, z_max;
            limitsForZ(m_Radius, &z_min, &z_max, &n, false);

            for(int j = 1; j < n; j++)
            {
                double z = z_min + j*(z_max -z_min)/(n - 1);
                double r = m_Radius;
                double dh(qMin(0.000625,r));

                int triangle_1 = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
                int triangle_2 = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z).shifted(-dh, 0));
                double temperature_1 = fem::triangles::functionValue(triangle_1, m_Net, m_T, geometry::_2d::Point(r, z));
                double temperature_2 = fem::triangles::functionValue(triangle_2, m_Net, m_T, geometry::_2d::Point(r, z).shifted(-dh, 0));
                //if(temperature_1 > 1685 + m_FrontHalfWidth && temperature_2 > 1685 + m_FrontHalfWidth)
                ostream << m_Materials[m_Net.material(triangle_1)].property(Material::HeatConductivity)*(temperature_1 - temperature_2)/dh << " " << z << endl;
            }
        }

        //************************************************************************************************************************
        //Распределение осевого потока тепла в радиальном направлении (прямая над фронтом).
        //************************************************************************************************************************
        if(Settings::instance()->setting("LocalConductiveHeatFlow/FlowAboveFront").toBool())
        {
            QFile file(qz_r_s_out_path.absolutePath() + "/" + filename);
            QTextStream ostream(&file);
            file.open(QFile::WriteOnly | QFile::Text);

            for(int i = 0; i < front_plus.size(); i++)
            {
                double dh(m_DiffStep);
                double r(front_plus[i].x());
                double z(front_plus[Settings::instance()->setting("LocalConductiveHeatFlow/CurvedSectionAboveFront").toBool()? i : 0].y() + 0.00125);

                int triangle_1 = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
                int triangle_2 = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z).shifted(0, dh));
                double temperature_1 = fem::triangles::functionValue(triangle_1, m_Net, m_T, geometry::_2d::Point(r, z));
                double temperature_2 = fem::triangles::functionValue(triangle_2, m_Net, m_T, geometry::_2d::Point(r, z).shifted(0, dh));

                if(m_Net.material(triangle_1) == m_Net.material(triangle_2))
                {
                    int material = m_Net.material(triangle_1);
                    ostream << r << " " << m_Materials[material].property(Material::HeatConductivity)*(temperature_2 - temperature_1)/dh << endl;
                }
            }
            file.close();
        }
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::temperatureProfiles()
{
    //****************************************************************************************************************************
    // Чтение и обработка настроек.
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_path(m_InputDir);
    // Маска файла/файлов для обработки.
    QStringList filename_mask(Settings::instance()->setting("TemperatureProfiles/FileMask").toString() + "T.dat");
    // Директория для записи результата.
    QDir out_path(m_OutputDir);

    // Получаем список файлов для обработки.
    QStringList filenames = in_path.entryList(filename_mask, QDir::CaseSensitive | QDir::Files);

    //****************************************************************************************************************************
    // Получение задаваемых извне координат сечений.
    //****************************************************************************************************************************
    QVariantList r_sections_temp = Settings::instance()->setting("TemperatureProfiles/R").toList();
    QVariantList z_sections_temp = Settings::instance()->setting("TemperatureProfiles/Z").toList();

    QList<double> r_sections;
    QList<double> z_sections;

    for(int i = 0; i < r_sections_temp.size(); i++)
        r_sections.append(r_sections_temp[i].toDouble());

    for(int i = 0; i < z_sections_temp.size(); i++)
        z_sections.append(z_sections_temp[i].toDouble());

    //****************************************************************************************************************************
    // Обработка файлов. Генерация профилей температуры.
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Температурные профили") << endl;
    foreach(QString filename, filenames)
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        m_T.readFromFile(in_path.absolutePath() + "/" + filename, linal::Vector::FEM_T);
        m_Net.readFromFiles(in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "rz.dat",
                            in_path.absolutePath() + "/" + filename.left(filename.size() - 5) + "nvtr.dat", fem::TriangleNet::NET_T);

        //****************************************************************************************************************************
        //Анализируем сетку и составляем таблицу, сокращающую поиск по сетке.
        //****************************************************************************************************************************
        m_NetSearcher->analyzeNet(m_Net);

        //****************************************************************************************************************************
        //Создаем профиль температуры в сечениях r=const.
        //****************************************************************************************************************************
        foreach(double r, r_sections)
        {
            if(!out_path.exists("T(z)/" + QString::number(r)))
                out_path.mkpath("T(z)/" + QString::number(r));

            QFile file(out_path.absolutePath() + "/T(z)/" + QString::number(r) + "/" + filename);
            QTextStream ostream(&file);

            int n;
            double z_min, z_max;
            limitsForZ(r, &z_min, &z_max, &n, false);

            file.open(QFile::WriteOnly | QFile::Text);
            for(int j = 0; j < n; j++)
            {
                double z = z_min + j*(z_max - z_min)/(n - 1);
                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
                ostream << fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(r, z)) << " " << z << "\n";
            }
            file.close();
        }

        //****************************************************************************************************************************
        //Создаем профиль температуры в сечениях z=const.
        //****************************************************************************************************************************
        foreach(double z, z_sections)
        {
            if(!out_path.exists("T(r)/" + QString::number(z)))
                out_path.mkpath("T(r)/" + QString::number(z));

            QFile file(out_path.absolutePath() + "/T(r)/" + QString::number(z) + "/" + filename);
            QTextStream ostream(&file);

            int n;
            double r_min, r_max;
            limitsForR(z, &r_min, &r_max, &n, false);

            file.open(QFile::WriteOnly | QFile::Text);
            for(int i = 0; i < n; i++)
            {
                double r = r_min + i*(r_max - r_min)/(n - 1);
                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
                ostream << r << " " << fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(r, z)) << "\n";
            }
            file.close();
        }
    }
    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::speedProfiles()
{
    //****************************************************************************************************************************
    // Чтение и обработка настроек.
    //****************************************************************************************************************************
    // Директория с входными данными.
    QDir in_path(m_InputDir);
    // Маска файла/файлов для обработки.
    QStringList filename_mask(Settings::instance()->setting("SpeedProfiles/FileMask").toString() + "V?.dat");
    // Директория для записи результата.
    QDir out_path(m_OutputDir);

    // Получаем список файлов для обработки.
    QStringList filenames = in_path.entryList(filename_mask, QDir::CaseSensitive | QDir::Files);

    //****************************************************************************************************************************
    // Получение задаваемых извне координат сечений.
    //****************************************************************************************************************************
    QVariantList r_sections_temp = Settings::instance()->setting("SpeedProfiles/R").toList();
    QVariantList z_sections_temp = Settings::instance()->setting("SpeedProfiles/Z").toList();

    QList<double> r_sections;
    QList<double> z_sections;

    for(int i = 0; i < r_sections_temp.size(); i++)
        r_sections.append(r_sections_temp[i].toDouble());

    for(int i = 0; i < z_sections_temp.size(); i++)
        z_sections.append(z_sections_temp[i].toDouble());

    //****************************************************************************************************************************
    // Обработка файлов. Генерация профилей скорости.
    //****************************************************************************************************************************
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    cout << QObject::tr("Профили скорости") << endl;

    foreach(QString filename, filenames)
    {
        cout << QObject::tr("Файл") << " " << filename << "\r";
        cout.flush();

        if(filename.contains("Vr"))
            profileVr(r_sections, filename, in_path.absolutePath(), out_path.absolutePath());
        if(filename.contains("Vz"))
            profileVz(z_sections, filename, in_path.absolutePath(), out_path.absolutePath());
    }

    cout << endl << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::profileVr(QList<double> r_sections, const QString &filename, const QString &in_path, const QString &out_path)
{
    m_Vr.readFromFile(in_path + "/" + filename, linal::Vector::FEM_T);
    m_Net.readFromFiles(in_path + "/" + filename.left(filename.size() - 6) + "rz.dat",
                        in_path + "/" + filename.left(filename.size() - 6) + "nvtr.dat", fem::TriangleNet::NET_T);

    //****************************************************************************************************************************
    //Анализируем сетку и составляем таблицу, сокращающую поиск по сетке.
    //****************************************************************************************************************************
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //Создаем профиль радиальной компоненты скорости.
    //****************************************************************************************************************************
    foreach(double r, r_sections)
    {
        if(!QDir(out_path).exists("Vr/" + QString::number(r)))
            QDir(out_path).mkpath("Vr/" + QString::number(r));

        QFile file(out_path + "/Vr/" + QString::number(r) + "/" + filename);
        QTextStream ostream(&file);

        int n;
        double z_min, z_max;
        limitsForZ(r, &z_min, &z_max, &n, true);

        file.open(QFile::WriteOnly | QFile::Text);
        for(int j = 0; j < n; j++)
        {
            double z = z_min + j*(z_max - z_min)/(n - 1);
            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            ostream << fem::triangles::functionValue(triangle, m_Net, m_Vr, geometry::_2d::Point(r, z)) << " " << z << endl;
        }
        file.close();
    }

    //****************************************************************************************************************************
    // Создаем профиль через центр конвективной ячейки.
    //****************************************************************************************************************************
    if(Settings::instance()->setting("SpeedProfiles/CentralSections").toBool())
    {
        m_Psi.readFromFile(in_path + "/" + filename.left(filename.size() - 6) + "Psi.dat", linal::Vector::FEM_T);

        geometry::_2d::Point centre = fem::triangles::minimum(m_Net, m_Psi);

        if(!QDir(out_path).exists("Vr/centre"))
            QDir(out_path).mkpath("Vr/centre");

        QFile file(out_path + "/Vr/centre/" + filename);
        QTextStream ostream(&file);

        int n;
        double z_min, z_max;
        limitsForZ(centre.x(), &z_min, &z_max, &n, true);

        file.open(QFile::WriteOnly | QFile::Text);
        for(int j = 0; j < n; j++)
        {
            double z = z_min + j*(z_max - z_min)/(n - 1);
            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(centre.x(), z));
            ostream << fem::triangles::functionValue(triangle, m_Net, m_Vr, geometry::_2d::Point(centre.x(), z)) << " " << z << endl;
        }
        file.close();
    }
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::profileVz(QList<double> z_sections, const QString &filename, const QString &in_path, const QString &out_path)
{
    m_Vz.readFromFile(in_path + "/" + filename, linal::Vector::FEM_T);
    m_Net.readFromFiles(in_path + "/" + filename.left(filename.size() - 6) + "rz.dat",
                        in_path + "/" + filename.left(filename.size() - 6) + "nvtr.dat", fem::TriangleNet::NET_T);

    //****************************************************************************************************************************
    //Анализируем сетку и составляем таблицу, сокращающую поиск по сетке.
    //****************************************************************************************************************************
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //Создаем профиль радиальной компоненты скорости.
    //****************************************************************************************************************************
    foreach(double z, z_sections)
    {
        if(!QDir(out_path).exists("Vz/" + QString::number(z)))
            QDir(out_path).mkpath("Vz/" + QString::number(z));

        QFile file(out_path + "/Vz/" + QString::number(z) + "/" + filename);
        QTextStream ostream(&file);

        int n;
        double r_min, r_max;
        limitsForR(z, &r_min, &r_max, &n, true);

        file.open(QFile::WriteOnly | QFile::Text);
        for(int i = 0; i < n; i++)
        {
            double r = r_min + i*(r_max - r_min)/(n - 1);
            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            ostream << r << " " << fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(r, z)) << endl;
        }
        file.close();
    }

    //************************************************************************************************************************
    //Распределение осевого потока тепла в радиальном направлении (прямая над фронтом).
    //************************************************************************************************************************
    if(Settings::instance()->setting("SpeedProfiles/VzAboveFront").toBool())
    {
        m_T.readFromFile(in_path + "/" + filename.left(filename.size() - 6) + "T.dat", linal::Vector::FEM_T);
        //********************************************************************************************************************
        // Выделяем фронт кристаллизации.
        //********************************************************************************************************************
        QList<geometry::_2d::Point> front_plus = identifyFrontExt(1685 + m_FrontHalfWidth);

        //********************************************************************************************************************
        // Выводим результат.
        //********************************************************************************************************************

        if(!QDir(out_path).exists("Vz/front"))
            QDir(out_path).mkpath("Vz/front");

        QFile file(out_path + "/Vz/front/" + filename);
        QTextStream ostream(&file);
        file.open(QFile::WriteOnly | QFile::Text);

        for(int i = 0; i < front_plus.size(); i++)
        {
            double r(front_plus[i].x());
            double z(front_plus[Settings::instance()->setting("SpeedProfiles/CurvedSectionAboveFront").toBool()? i : 0].y() + 0.00125);

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            ostream << r << " " << fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(r, z)) << endl;
        }
        file.close();
    }

    //****************************************************************************************************************************
    // Создаем профиль через центр конвективной ячейки.
    //****************************************************************************************************************************
    if(Settings::instance()->setting("SpeedProfiles/CentralSections").toBool())
    {
        m_Psi.readFromFile(in_path + "/" + filename.left(filename.size() - 6) + "Psi.dat", linal::Vector::FEM_T);

        geometry::_2d::Point centre = fem::triangles::minimum(m_Net, m_Psi);

        if(!QDir(out_path).exists("Vz/centre"))
            QDir(out_path).mkpath("Vz/centre");

        QFile file(out_path + "/Vz/centre/" + filename);
        QTextStream ostream(&file);
        file.open(QFile::WriteOnly | QFile::Text);

        int n;
        double r_min, r_max;
        limitsForR(centre.y(), &r_min, &r_max, &n, true);

        for(int i = 0; i < n; i++)
        {
            double r = r_min + i*(r_max - r_min)/(n - 1);
            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, centre.y()));
            ostream << r << " " << fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(r, centre.y())) << endl;
        }
        file.close();

    }
}

//================================================================================================================================
//================================================================================================================================
bool AbstractSolver::solveSLAE(const linal::CSRMatrix &A, const linal::Vector &f, linal::Vector &x)
{
    foreach(linal::SLAESolver* solver, m_Solvers)
        if(solver->solveLU(A,f,x) < m_EpsSLAE)
            return true;

    foreach(linal::SLAESolver* solver, m_Solvers)
        if(solver->solveDi(A,f,x) < m_EpsSLAE)
            return true;

    foreach(linal::SLAESolver* solver, m_Solvers)
        if(solver->solve(A,f,x) < m_EpsSLAE)
            return true;

    qWarning("SLAE wasn't solved!");
    return false;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::solve()
{
    if(m_SubTasks.testFlag(AbstractSolver::Stationarity))
        solveStationary();
    else
        solveNonStationary();
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::solveStationary()
{
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    if(!m_SubTasks.testFlag(AbstractSolver::Convection))
    {
        //Решаем уравнение для T
        m_A.zeroize();
        m_F.zeroize();
        formTSLAE(m_TimeMoment);
        solveSLAE(m_A, m_F, m_T);
        writeResult(m_TimeMoment);
    }

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    else if(m_SubTasks.testFlag(AbstractSolver::Convection))
    {
        int nolin_iter(0);
        double nolin_eps(0.);

        do
        {
            formV();

            //Решаем уравнение для T
            m_A.zeroize();
            m_F.zeroize();
            formTSLAE(m_TimeMoment);
            solveSLAE(m_A, m_F, m_T);
            m_T = (1 - m_W)*m_TOldNolin + m_W*m_T;

            //Решаем уравнение для Vf
            if(m_SubTasks.testFlag(AbstractSolver::Rotation))
            {
                m_A.zeroize();
                m_F.zeroize();
                formVfSLAE(m_TimeMoment);
                solveSLAE(m_A, m_F, m_Vf);
                m_Vf = (1 - m_W)*m_VfOldNolin + m_W*m_Vf;
            }

            //Решаем уравнение для Omega
            m_A.zeroize();
            m_F.zeroize();
            formOmegaSLAE(m_TimeMoment);
            solveSLAE(m_A, m_F, m_Omega);
            m_Omega = (1 - m_W)*m_OmegaOldNolin + m_W*m_Omega;

            //Решаем уравнение для Psi
            m_A.zeroize();
            m_F.zeroize();
            formPsiSLAE(m_TimeMoment);
            solveSLAE(m_A, m_F, m_Psi);
            m_Psi = (1 - m_W)*m_PsiOldNolin + m_W*m_Psi;

            //Рассчитываем невязку по нелинейности
            double res_t =linal::Vector::normSquared(m_T - m_TOldNolin)/linal::Vector::normSquared(m_T);
            double res_psi = linal::Vector::normSquared(m_Psi - m_PsiOldNolin)/linal::Vector::normSquared(m_Psi);
            double res_omega = linal::Vector::normSquared(m_Omega - m_OmegaOldNolin)/linal::Vector::normSquared(m_Omega);
            nolin_eps = qMax(qSqrt(res_t)/m_W, qMax(qSqrt(res_psi)/m_W, qSqrt(res_omega)/m_W));

            if(m_SubTasks.testFlag(AbstractSolver::Rotation))
            {
                double res_vf = linal::Vector::normSquared(m_Vf - m_VfOldNolin)/linal::Vector::normSquared(m_Vf);
                nolin_eps = qMax(nolin_eps, qSqrt(res_vf)/m_W);

                m_VfOldNolin = m_Vf;
            }

            m_TOldNolin = m_T;
            m_PsiOldNolin = m_Psi;
            m_OmegaOldNolin = m_Omega;
        }
        while(++nolin_iter < m_MaxiterNonlin && nolin_eps > m_EpsNonlin);
        writeResult(m_TimeMoment);
    }
    else
        cout << QObject::tr("Неизвестная задача") << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::solveNonStationary()
{
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    if(m_ResumeFrom != -1)
        readResumeInfo();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    if(!m_SubTasks.testFlag(AbstractSolver::Convection))
    {
        int first_layer = (m_ResumeFrom == -1 ? 1 : m_ResumeFrom + 1);
        for(int cur_layer = first_layer; cur_layer < m_NumOfTimeLayers; cur_layer++)
        {
            //Решаем уравнение для T
            m_A.zeroize();
            m_F.zeroize();
            formTSLAE(cur_layer);
            solveSLAE(m_A, m_F, m_T);

            if(m_SubTasks.testFlag(AbstractSolver::Solidification))
                recreateNet();

            if(m_SubTasks.testFlag(AbstractSolver::Concentration))
            {
                m_A.zeroize();
                m_F.zeroize();
                formCSLAE(cur_layer);
                solveSLAE(m_A, m_F, m_C);

                m_COldTime = m_C;

                QList<double> args, vals, points, rez;
                for(int i = 0; i < m_FrontPlus.size(); i++)
                {
                    args.append(m_Net.point(m_FrontPlus[i]).x());
                    vals.append(m_C[m_FrontPlus[i]]);
                }
                for(int i = 0; i < m_FrontMinus.size(); i++)
                    points.append(m_Net.point(m_FrontMinus[i]).x());

                rez = approximation::Interpolation::interpolate(args, vals, points, approximation::Interpolation::Linear);
                for(int i = 0; i < m_FrontMinus.size(); i++)
                    m_C[m_FrontMinus[i]] = rez[i]*m_Materials[5].property(Material::SegregationFactor);
            }

            m_TOldTime = m_T;

            if(cur_layer%m_OutputStep == 0)
            {
                writeResult(cur_layer);
                if(m_SubTasks.testFlag(AbstractSolver::Debug))
                    writeResumeInfo(cur_layer);
            }
            if(cur_layer%50 == 0)
                cout << cur_layer << endl;
        }
    }
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    else if(m_SubTasks.testFlag(AbstractSolver::Convection))
    {
        int first_layer = (m_ResumeFrom == -1 ? 1 : m_ResumeFrom + 1);
        for(int cur_layer = first_layer; cur_layer < m_NumOfTimeLayers; cur_layer++)
        {
            int nolin_iter(0);
            double nolin_eps(0.);

            do
            {
                formV();

                //Решаем уравнение для T
                m_A.zeroize();
                m_F.zeroize();
                formTSLAE(cur_layer);
                solveSLAE(m_A, m_F, m_T);
                m_T = (1 - m_W)*m_TOldNolin + m_W*m_T;

                if(m_SubTasks.testFlag(AbstractSolver::Solidification))
                    recreateNet();
                //writeResult(cur_layer);

                //Решем уравнения для Vf
                if(m_SubTasks.testFlag(AbstractSolver::Rotation))
                {
                    m_A.zeroize();
                    m_F.zeroize();
                    formVfSLAE(cur_layer);
                    solveSLAE(m_A, m_F, m_Vf);
                    m_Vf = (1 - m_W)*m_VfOldNolin + m_W*m_Vf;
                }

                //Решаем уравнение для Omega
                m_A.zeroize();
                m_F.zeroize();
                formOmegaSLAE(cur_layer);
                solveSLAE(m_A, m_F, m_Omega);
                m_Omega = (1 - m_W)*m_OmegaOldNolin + m_W*m_Omega;

                //Решаем уравнение для Psi
                m_A.zeroize();
                m_F.zeroize();
                formPsiSLAE(cur_layer);
                solveSLAE(m_A, m_F, m_Psi);
                m_Psi = (1 - m_W)*m_PsiOldNolin + m_W*m_Psi;

                //Рассчитываем невязку по нелинейности
                double res_t = linal::Vector::normSquared(m_T - m_TOldNolin)/linal::Vector::normSquared(m_T);
                double res_psi = linal::Vector::normSquared(m_Psi - m_PsiOldNolin)/linal::Vector::normSquared(m_Psi);
                double res_omega = linal::Vector::normSquared(m_Omega - m_OmegaOldNolin)/linal::Vector::normSquared(m_Omega);
                nolin_eps = qMax(qSqrt(res_t)/m_W, qMax(qSqrt(res_psi)/m_W, qSqrt(res_omega)/m_W));

                if(m_SubTasks.testFlag(AbstractSolver::Rotation))
                {
                    double res_vf = linal::Vector::normSquared(m_Vf - m_VfOldNolin)/linal::Vector::normSquared(m_Vf);
                    nolin_eps = qMax(nolin_eps, qSqrt(res_vf)/m_W);

                    m_VfOldNolin = m_Vf;
                }

                m_TOldNolin = m_T;
                m_PsiOldNolin = m_Psi;
                m_OmegaOldNolin = m_Omega;
            }
            while(++nolin_iter < m_MaxiterNonlin && nolin_eps > m_EpsNonlin);

            if(m_SubTasks.testFlag(AbstractSolver::Concentration))
            {
                m_A.zeroize();
                m_F.zeroize();
                formCSLAE(cur_layer);
                solveSLAE(m_A, m_F, m_T);

                m_COldTime = m_C;

                QList<double> args, vals, points, rez;
                for(int i = 0; i < m_FrontPlus.size(); i++)
                {
                    args.append(m_Net.point(m_FrontPlus[i]).x());
                    vals.append(m_C[m_FrontPlus[i]]);
                }
                for(int i = 0; i < m_FrontMinus.size(); i++)
                    points.append(m_Net.point(m_FrontMinus[i]).x());

                rez = approximation::Interpolation::interpolate(args, vals, points, approximation::Interpolation::Linear);
                for(int i = 0; i < m_FrontMinus.size(); i++)
                    m_C[m_FrontMinus[i]] = rez[i]*m_Materials[5].property(Material::SegregationFactor);
            }

            m_TOldTime = m_T;
            m_OmegaOldTime = m_Omega;
            if(m_SubTasks.testFlag(AbstractSolver::Rotation))
                m_VfOldTime = m_Vf;

            if(cur_layer%m_OutputStep == 0)
            {
                writeResult(cur_layer);
                if(m_SubTasks.testFlag(AbstractSolver::Debug))
                    writeResumeInfo(cur_layer);
            }

            cout << cur_layer << ": " << nolin_iter << " " << nolin_eps << endl;
        }
    }
    else
        cout << QObject::tr("Неизвестная задача") << endl;
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::formTSLAE(int time_layer)
{
    double density;			// Нижние подчеркивания для того чтобы имя переменной не совпадало с именем функции.
    double heat_capacity;
    double heat_conductivity;
    double D[3][3];     		//Матрица, используемая для определения базисных функций
    double C[3][3][3];			//Локальная матрица массы
    double detD;        		//Определитель матрицы D
    double abs_detD;    		//Модуль определетиля матрицы D
    double ht;          		//Шаг по времени
    double r[3];        		//r-координаты вершин КЭ
    double z[3];        		//z-координаты вершин КЭ
    int top[3];         		//Индексы вершин КЭ

    ht = m_TimeStep;			// Задел на будущее, когда сетка может быть неравномерной.

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        density = m_Materials[m_Net.material(cur_triangle)].property(Material::Density);
        heat_capacity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatCapacity);
        heat_conductivity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatConductivity);

        for(int i = 0; i < 3; i++)
        {
            r[i] = m_Net.point(cur_triangle, i).x();
            z[i] = m_Net.point(cur_triangle, i).y();
            top[i] = m_Net.pointIndex(cur_triangle, i);
        }

        detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
        abs_detD = qAbs(detD);

        D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
        D[0][1] = (z[1] - z[2])/detD;
        D[0][2] = (r[2] - r[1])/detD;
        D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
        D[1][1] = (z[2] - z[0])/detD;
        D[1][2] = (r[0] - r[2])/detD;
        D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
        D[2][1] = (z[0] - z[1])/detD;
        D[2][2] = (r[1] - r[0])/detD;

        C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
        C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
        C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

        C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
        C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
        C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

        C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
        C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
        C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                double addition(0);
                addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*heat_conductivity*abs_detD/6.;
                if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                    addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])*heat_capacity*density/ht;
                if(m_SubTasks.testFlag(AbstractSolver::Convection))
                {
                    addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1]*heat_capacity*density;
                    addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2]*heat_capacity*density;
                }
                m_A(top[i], top[j]) += addition;
            }
        }

        if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        m_F[top[i]] += m_TOldTime[top[j]]*r[k]*C[i][j][k]*heat_capacity*density/ht;
    }

    //****************************************************************************************************************************
    //Учитываем краевые условия.
    //****************************************************************************************************************************
    tBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::formCSLAE(int time_layer)
{
    double diffusivity;
    double D[3][3];     		//Матрица, используемая для определения базисных функций
    double C[3][3][3];			//Локальная матрица массы
    double detD;        		//Определитель матрицы D
    double abs_detD;    		//Модуль определетиля матрицы D
    double ht;          		//Шаг по времени
    double r[3];        		//r-координаты вершин КЭ
    double z[3];        		//z-координаты вершин КЭ
    int top[3];         		//Индексы вершин КЭ

    ht = m_TimeStep;			// Задел на будущее, когда сетка может быть неравномерной.

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            diffusivity = m_Materials[5].property(Material::Diffusivity);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*diffusivity*abs_detD/6.;
                    if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                        addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])/ht;
                    if(m_SubTasks.testFlag(AbstractSolver::Convection))
                    {
                        addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                     m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                     m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1];
                        addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                     m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                     m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2];
                    }
                    m_A(top[i], top[j]) += addition;
                }
            }

            if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                for(int i = 0; i < 3; i++)
                    for(int j = 0; j < 3; j++)
                        for(int k = 0; k < 3; k++)
                            m_F[top[i]] += m_COldTime[top[j]]*r[k]*C[i][j][k]/ht;
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    cBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::formPsiSLAE(int time_layer)
{
    double D[3][3];     //Матрица, используемая для определения базисных функций
    double C[3][3][3];	//Локальная матрица массы
    double detD;        //Определитель матрицы D
    double abs_detD;    //Модуль определетиля матрицы D
    double r[3];        //r-координаты вершин КЭ
    double z[3];        //z-координаты вершин КЭ
    int top[3];         //Индексы вершин КЭ

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_A(top[i], top[j]) += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*abs_detD/6.;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_F[top[i]] += m_Omega[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2]);

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_PsiOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])/r[j];
                    else
                    {
                        double psi_diff = derivation(m_PsiOldNolin, top[j], AbstractSolver::R, AbstractSolver::RDiff);
                        m_F[top[i]] -= psi_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2]);
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    psiBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::formVfSLAE(int time_layer)
{
    double viscosity;
    double D[3][3];     	//Матрица, используемая для определения базисных функций
    double C[3][3][3];		//Локальная матрица массы
    double detD;        	//Определитель матрицы D
    double abs_detD;    	//Модуль определетиля матрицы D
    double ht;          	//Шаг по времени
    double r[3];        	//r-координаты вершин КЭ
    double z[3];        	//z-координаты вершин КЭ
    int top[3];         	//Индексы вершин КЭ

    ht = m_TimeStep;		// Задел на будущее, когда сетка может быть неравномерной.

    //*******************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //*******************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;

            m_F[m_Net.pointIndex(cur_triangle, 0)] = m_Net.point(cur_triangle, 0).x()*m_VRot;
            m_F[m_Net.pointIndex(cur_triangle, 1)] = m_Net.point(cur_triangle, 1).x()*m_VRot;
            m_F[m_Net.pointIndex(cur_triangle, 2)] = m_Net.point(cur_triangle, 2).x()*m_VRot;
        }
        else
        {
            viscosity = m_Materials[m_Net.material(cur_triangle)].property(Material::Viscosity);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*viscosity*abs_detD/6.;
                    if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                        addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])/ht;
                    addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1];
                    addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2];
                    addition += m_Vr[top[0]]*C[i][j][0];
                    addition += m_Vr[top[1]]*C[i][j][1];
                    addition += m_Vr[top[2]]*C[i][j][2];

                    m_A(top[i], top[j]) += addition;
                }
            }

            if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                for(int i = 0; i < 3; i++)
                    for(int j = 0; j < 3; j++)
                        m_F[top[i]] += m_VfOldTime[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2])/ht;


            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_VfOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity/r[j];
                    else
                    {
                        double omega_diff = derivation(m_VfOldNolin, top[j], AbstractSolver::R, AbstractSolver::RDiff);
                        m_F[top[i]] -= omega_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity;
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    vfBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void AbstractSolver::formOmegaSLAE(int time_layer)
{
    double viscosity;		// Нижние подчеркивания для того чтобы имя переменной не совпадало с именем функции.
    double expansiveness;	// Нижние подчеркивания для того чтобы имя переменной не совпадало с именем функции.
    double D[3][3];     	//Матрица, используемая для определения базисных функций
    double C[3][3][3];		//Локальная матрица массы
    double detD;        	//Определитель матрицы D
    double abs_detD;    	//Модуль определетиля матрицы D
    double ht;          	//Шаг по времени
    double r[3];        	//r-координаты вершин КЭ
    double z[3];        	//z-координаты вершин КЭ
    int top[3];         	//Индексы вершин КЭ

    ht = m_TimeStep;		// Задел на будущее, когда сетка может быть неравномерной.

    //*******************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //*******************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            viscosity = m_Materials[m_Net.material(cur_triangle)].property(Material::Viscosity);
            expansiveness = m_Materials[m_Net.material(cur_triangle)].property(Material::Expansiveness);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*viscosity*abs_detD/6.;
                    if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                        addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])/ht;
                    addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1];
                    addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2];
                    addition -= m_Vr[top[0]]*C[i][j][0];
                    addition -= m_Vr[top[1]]*C[i][j][1];
                    addition -= m_Vr[top[2]]*C[i][j][2];

                    m_A(top[i], top[j]) += addition;
                }
            }

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        m_F[top[i]] -= m_T[top[j]]*r[k]*(C[i][k][0] + C[i][k][1] + C[i][k][2])*D[j][1]*expansiveness*9.8;

            if(!m_SubTasks.testFlag(AbstractSolver::Stationarity))
                for(int i = 0; i < 3; i++)
                    for(int j = 0; j < 3; j++)
                        m_F[top[i]] += m_OmegaOldTime[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2])/ht;

            if(m_SubTasks.testFlag(AbstractSolver::Rotation))
                for(int i = 0; i < 3; i++)
                    for(int j = 0; j < 3; j++)
                        for(int k = 0; k < 3; k++)
                            m_F[top[i]] += 2*m_Vf[top[j]]*m_Vf[top[k]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*D[k][2];

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_OmegaOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity/r[j];
                    else
                    {
                        double omega_diff = derivation(m_OmegaOldNolin, top[j],  AbstractSolver::R, AbstractSolver::RDiff);
                        m_F[top[i]] -= omega_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity;
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    omegaBC(time_layer);
}

//================================================================================================================================
//Вывод результата в файл.
//================================================================================================================================
void AbstractSolver::writeResult(int time_layer)
{
    QFile file;
    QDataStream ostream(&file);
    ostream.setVersion(QDataStream::Qt_4_7);

    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    //****************************************************************************************************************************
    //Выводим конечноэлементую сетку.
    //****************************************************************************************************************************
    QString filename_nvtr = QString("/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_rz = QString("/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    m_Net.writeToFiles(filename_rz, filename_nvtr);

    //****************************************************************************************************************************
    //Выводим значения температуры.
    //****************************************************************************************************************************
    file.setFileName(QString("/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
    file.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_T[i];
    file.close();

    if(m_SubTasks.testFlag(AbstractSolver::Concentration))
    {
        //****************************************************************************************************************************
        //Выводим значения концентрации.
        //****************************************************************************************************************************
        file.setFileName(QString("/%1C.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
        file.open(QFile::WriteOnly);
        for(int i = 0; i < m_Net.numOfPoints(); i++)
            ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_C[i];
        file.close();
    }

    if(m_SubTasks.testFlag(AbstractSolver::Convection))
    {
        //************************************************************************************************************************
        //Выводим значения функции тока.
        //************************************************************************************************************************
        file.setFileName(QString("/%1Psi.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
        file.open(QFile::WriteOnly);
        for(int i = 0; i < m_Net.numOfPoints(); i++)
            ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_Psi[i];
        file.close();

        //************************************************************************************************************************
        //Выводим значения вихря.
        //************************************************************************************************************************
        file.setFileName(QString("/%1Omega.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
        file.open(QFile::WriteOnly);
        for(int i = 0; i < m_Net.numOfPoints(); i++)
            ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_Omega[i];
        file.close();

        //************************************************************************************************************************
        //Выводим значения r-компоненты скорости.
        //************************************************************************************************************************
        file.setFileName(QString("/%1Vr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
        file.open(QFile::WriteOnly);
        for(int i = 0; i < m_Net.numOfPoints(); i++)
            ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_Vr[i];
        file.close();

        //************************************************************************************************************************
        //выводим значения z-компоненты скорости.
        //************************************************************************************************************************
        file.setFileName(QString("/%1Vz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
        file.open(QFile::WriteOnly);
        for(int i = 0; i < m_Net.numOfPoints(); i++)
            ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_Vz[i];
        file.close();

        if(m_SubTasks.testFlag(AbstractSolver::Rotation))
        {
            //********************************************************************************************************************
            //выводим значения f-компоненты скорости.
            //********************************************************************************************************************
            file.setFileName(QString("/%1Vf.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputDir));
            file.open(QFile::WriteOnly);
            for(int i = 0; i < m_Net.numOfPoints(); i++)
                ostream << m_Net.point(i).x() << m_Net.point(i).y() << m_Vf[i];
            file.close();
        }
    }
}

//================================================================================================================================
//Инициализация структур данных, необходимая для продолжения расчетов с указанного места.
//================================================================================================================================
void AbstractSolver::readResumeInfo()
{
    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    QString filename_rz = QString("/%1rz.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_nvtr = QString("/%1nvtr.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);

    QString filename_t = QString("/%1T.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_c = QString("/%1C.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_vf = QString("/%1Vf.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_psi = QString("/%1Psi.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_omega = QString("/%1Omega.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);

    QString filename_bc = QString("/%1bc.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);
    QString filename_net = QString("/%1net.dat").arg(m_ResumeFrom, num_of_orders, 10, QChar('0')).prepend(m_OutputDir);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    if(Settings::instance()->setting("BinaryInput").toBool())
        m_Net.readFromFiles(filename_rz, filename_nvtr, fem::TriangleNet::NET_B);
    else
        m_Net.readFromFiles(filename_rz, filename_nvtr, fem::TriangleNet::NET_T);
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int data_size;
    int bc_top, bc_type;

    QFile file(filename_bc);
    QTextStream istream(&file);
    file.open(QIODevice::ReadOnly | QIODevice::Text);

    istream >> data_size;
    m_TopTBC.clear();
    m_TypeTBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopTBC.append(bc_top);
        m_TypeTBC.append(bc_type);
    }

    istream >> data_size;
    m_TopPsiBC.clear();
    m_TypePsiBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopPsiBC.append(bc_top);
        m_TypePsiBC.append(bc_type);
    }

    istream >> data_size;
    m_TopOmegaBC.clear();
    m_TypeOmegaBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopOmegaBC.append(bc_top);
        m_TypeOmegaBC.append(bc_type);
    }

    istream >> data_size;
    m_TopVfBC.clear();
    m_TypeVfBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopVfBC.append(bc_top);
        m_TypeVfBC.append(bc_type);
    }
    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    file.setFileName(filename_net);
    file.open(QIODevice::ReadOnly | QIODevice::Text);
    switch(istream.readLine().toInt())
    {
    case 0: m_State = AbstractSolver::NotStarted; break;
    case 1: m_State = AbstractSolver::JustStarted; break;
    case 2: m_State = AbstractSolver::InProgress; break;
    case 3: m_State = AbstractSolver::Finished; break;
    default: qFatal("Unknown crystallization state");
    }
    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_T.resize(m_Net.numOfPoints());
    m_TOldTime.resize(m_Net.numOfPoints());
    m_TOldNolin.resize(m_Net.numOfPoints());

    linal::Vector::FileFormat format;
    if(Settings::instance()->setting("BinaryInput").toBool())
        format = linal::Vector::FEM_B;
    else
        format = linal::Vector::FEM_T;

    m_T.readFromFile(filename_t, format);
    m_TOldTime = m_T;
    m_TOldNolin = m_T;

    if(m_SubTasks.testFlag(AbstractSolver::Concentration) == true)
    {
        m_C.readFromFile(filename_c, format);
        m_COldTime = m_C;
    }

    if(m_SubTasks.testFlag(AbstractSolver::Convection) == true)
    {
        m_Psi.readFromFile(filename_psi, format);
        m_PsiOldNolin = m_Psi;

        m_Omega.readFromFile(filename_omega, format);
        m_OmegaOldTime = m_Omega;
        m_OmegaOldNolin = m_Omega;

        if(m_SubTasks.testFlag(AbstractSolver::Rotation))
        {
            m_Vf.readFromFile(filename_vf, format);
            m_VfOldTime = m_Vf;
            m_VfOldNolin = m_Vf;
        }
    }

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_F.resize(m_Net.numOfPoints());

    QVector<int> ig, jg, ig2, jg2;
    fem::triangles::formMatrixProfile(m_Net, ig, jg, ig2, jg2);
    m_A = linal::CSRMatrix(m_Net.numOfPoints(), ig, jg, ig2, jg2);
}

//================================================================================================================================
//Запись различной информации о расчетной области и краевых условиях нужной для возобновления расчетов с любого места.
//================================================================================================================================
void AbstractSolver::writeResumeInfo(int time_layer)
{
    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    QFile file(m_OutputDir % QString("/%1bc.dat").arg(time_layer, num_of_orders, 10, QChar('0')));
    QTextStream ostream(&file);

    file.open(QIODevice::WriteOnly | QIODevice::Text);
    ostream << m_TopTBC.size() << endl;
    for(int i = 0; i < m_TopTBC.size(); i++)
        ostream << m_TopTBC[i] << " " << m_TypeTBC[i] << endl;

    ostream << m_TopPsiBC.size() << endl;
    for(int i = 0; i < m_TopPsiBC.size(); i++)
        ostream << m_TopPsiBC[i] << " " << m_TypePsiBC[i] << endl;

    ostream << m_TopOmegaBC.size() << endl;
    for(int i = 0; i < m_TopOmegaBC.size(); i++)
        ostream << m_TopOmegaBC[i] << " " << m_TypeOmegaBC[i] << endl;

    ostream << m_TopVfBC.size() << endl;
    for(int i = 0; i < m_TopVfBC.size(); i++)
        ostream << m_TopVfBC[i] << " " << m_TypeVfBC[i] << endl;
    file.close();

    file.setFileName(m_OutputDir % QString("/%1net.dat").arg(time_layer, num_of_orders, 10, QChar('0')));
    file.open(QIODevice::WriteOnly | QIODevice::Text);
    ostream << m_State << endl;
    file.close();
}
